% 求k
clear,clc;
syms t x(t) v(t) k f(t) m g r;
cd1 = [m,g,r;300,9.8,3];%求k时用到的常数
f(t)=v*k*r^2;
eqn = [m*diff(v,t) == m*g-f];
cond = [v(0)==0];
[vt]=dsolve(eqn,cond);
xt = -1*int(vt,t,[0 t])+500; 
x_t_01 = subs(xt,cd1(1,:),cd1(2,:));
x_t_01 = matlabFunction(x_t_01);%得到表二下的x(t)
xdata = [0:3:30];
ydata = [500 470 425 372 317 264 215 160 108 55 1];
[resk,res] = lsqcurvefit(x_t_01,rand(1),xdata,ydata);

%测试k值
% for i = 1:length(xdata)
%     testk = x_t_01(resk,xdata(i));
%     vpa(testk)
% end
save fun xt vt m g r t resk k; 
clear;
load fun;
syms n;
cd2 = [m k g;2000/n resk 9.8];%求解r、n时用到的常数
global funcX funcV des;
funcX = subs(xt,cd2(1,:),cd2(2,:));
funcV = subs(vt,cd2(1,:),cd2(2,:));
des = [n r t];%决策变量

%拟合得到半径、价格曲线
r_choice = [2:0.5:4];
price =[65 170 350 660 1100]; 
[ab]=lsqcurvefit(@price_func,rand(2,1),r_choice,price);
test_x = [2:0.01:4];
test_y = ab(1)*test_x.^ab(2);
plot(test_x,test_y,'-',r_choice,price,'*');

%使用nlp。
x0 = rand(3,1); 
cost = @(x) x(1)*(ab(1)*x(2)^ab(2)+16*4*x(2)*sqrt(2)+200);
lb = [1 2 0];
ub = [100 4 100];

%使用非线性、整数混合规划。r取（4~8）/2
opti_cost = @(x) x(1)*(ab(1)*(x(2)/2)^ab(2)+16*4*(x(2)/2)*sqrt(2)+200);
%opti_cost = @(x) finalCost([x(1) x(2)/2 x(3)]);
opti_lb = [1 4 0];
opti_ub = [100 8 100];
fmin_x = fmincon(cost,x0,[],[],[],[],lb,ub,@nlc);
cl = [-inf 0];
cu = [0 0];
xtype = 'IIC';
opts = optiset('solver','BONMIN','display','iter');
Opt = opti('fun',opti_cost,'nl',@opti_nlc,cl,cu,'bounds',opti_lb,opti_ub,'xtype',xtype,'options',opts,'x0',x0);
[opti_x,~,ef,info]=solve(Opt,x0);
opti_x(2) = opti_x(2)/2;
fval = finalCost(opti_x);

function [g,h]=nlc(x)   
global funcX funcV des;
g(1) = double(subs(funcV,des,[x(1),x(2),x(3)])-20);
g(2) = double(subs(funcX,des,[x(1),x(2),x(3)])); 
h = [];
end

function [nlc]=opti_nlc(x)   
global funcX funcV des;
g = double(subs(funcV,des,[x(1),x(2)/2,x(3)])-20);
h = double(subs(funcX,des,[x(1),x(2)/2,x(3)])); 
nlc = [g h];
end

%用函数C1=a*r^bl来拟合r与单个伞面价格曲线
function [ydata] = price_func(x,xdata)
    ydata = x(1)*xdata.^x(2);
end

function [C] = finalCost(x)
C1 = 10000;%单个降落伞伞面费用
if x(2)<=2
    C1 = 65;
elseif 2<x(2)&&x(2)<=2.5
    C1 = 170;
elseif 2.5<x(2)&&x(2)<=3
    C1 = 350;
elseif  3<x(2)&&x(2)<=3.5
    C1 = 660;
elseif  3.5<x(2)&&x(2)<=4
    C1 = 1000;
end
%伞个数*(单个伞面价格+绳子数*单位长绳子价格*单条绳子长+单个伞固定费用)
C = x(1)*(C1+16*4*x(2)*sqrt(2)+200);
% dis = [x';C rg rh]
end











    
    
    
    



